function [u, w] = comp_velocity(c_top, h, p, sp)

p.Ma = 1;

k = sp.k;

c_hat = fft(c_top);

z = h * linspace(0, 1, sp.Nz);

for j = 1:sp.Nz
    w_hat = -0.4e1 * c_hat .* p.Ma .* ((((-0.1e1 / 0.4e1 + k * z(j) / 0.2e1) .* exp(0.2e1 * k * z(j)) + k .* (z(j) * h .* k + h / 0.2e1 - z(j) / 0.2e1)) .* exp(0.2e1 * k * h) + (k * h + 0.1e1 / 0.2e1) .* (-0.1e1 / 0.2e1 + k .* z(j)) .* exp(0.2e1 * k * z(j)) + k .* (h - z(j)) / 0.2e1) .* exp(-k * z(j)) - exp(k * z(j)) .* (-0.1e1 / 0.2e1 + (-0.1e1 / 0.2e1 + k * h) .* exp(0.2e1 * k * h)) / 0.2e1) .* exp(k * h) ./ (0.4e1 .* exp(0.4e1 * k * h) + 0.8e1 * exp(0.2e1 * k * h) + 0.16e2 * h ^ 2 * k .^ 2 .* exp(0.2e1 * k * h) + 0.4e1);
    w(j,:) = real(ifft(w_hat));
    u_hat = (-4*1i) * p.Ma * c_hat .* exp(k * h) .* ((((-0.1e1 / 0.4e1 + k * z(j) / 0.2e1) .* exp(0.2e1 * k * z(j)) - 0.1e1 / 0.2e1 - k .^ 2 * z(j) * h + (h / 0.2e1 + z(j) / 0.2e1) * k) .* exp(0.2e1 * k * h) + (k * h + 0.1e1 / 0.2e1) .* (-0.1e1 / 0.2e1 + k * z(j)) .* exp(0.2e1 * k * z(j)) - 0.1e1 / 0.2e1 + (-h / 0.2e1 + z(j) / 0.2e1) * k) .* exp(-k * z(j)) - exp(k * z(j)) .* ((k * h - 0.3e1 / 0.2e1) .* exp(0.2e1 * k * h) - 0.2e1 * k * h - 0.3e1 / 0.2e1) / 0.2e1) ./ (0.4e1 * exp(0.4e1 * k * h) + 0.8e1 * exp(0.2e1 * k * h) + 0.16e2 * h ^ 2 * k .^ 2 .* exp(0.2e1 * k * h) + 0.4e1);
    u(j,:) = real(ifft(u_hat));
end